home *** CD-ROM | disk | FTP | other *** search
/ System Booster / System Booster.iso / Archives / GNU / GNUPLOTsrc.lha / fit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-01-22  |  36.5 KB  |  1,385 lines

  1. #ifndef lint
  2. static char *RCSid = "$Id: fit.c,v 1.33 1995/12/11 23:14:04 drd Exp $";
  3. #endif
  4.  
  5. /*
  6.  *    Nonlinear least squares fit according to the
  7.  *    Marquardt-Levenberg-algorithm
  8.  *
  9.  *    added as Patch to Gnuplot (v3.2 and higher)
  10.  *    by Carsten Grammes
  11.  *    Experimental Physics, University of Saarbruecken, Germany
  12.  *
  13.  *    Internet address: cagr@rz.uni-sb.de
  14.  *
  15.  *    Copyright of this module:   Carsten Grammes, 1993
  16.  *
  17.  *    Permission to use, copy, and distribute this software and its
  18.  *    documentation for any purpose with or without fee is hereby granted,
  19.  *    provided that the above copyright notice appear in all copies and
  20.  *    that both that copyright notice and this permission notice appear
  21.  *    in supporting documentation.
  22.  *
  23.  *    This software is provided "as is" without express or implied warranty.
  24.  *
  25.  *    930726:     Recoding of the Unix-like raw console I/O routines by:
  26.  *            Michele Marziani (marziani@ferrara.infn.it)
  27.  * drd: start unitialised variables at 1 rather than NEARLY_ZERO
  28.  *  (fit is more likely to converge if started from 1 than 1e-30 ?)
  29.  */
  30.  
  31.  
  32. #define FIT_MAIN
  33.  
  34. #include <math.h>
  35. #include <ctype.h>
  36. #include <signal.h>
  37.  
  38. #if defined(MSDOS) || defined(DOS386)         /* non-blocking IO stuff */
  39. #include <io.h>
  40. #include <conio.h>
  41. #include <dos.h>
  42. #else
  43. #if defined(AMIGA_SC_6_1)
  44. #include <fcntl.h>
  45. #include <ios1.h>
  46. #include <exec/types.h>
  47. #include <exec/tasks.h>
  48. #include <dos/dos.h>
  49. #include <proto/exec.h>
  50. #include <proto/dos.h>
  51. #else
  52. #ifdef OSK
  53. #include <stdio.h>
  54. #else
  55. #ifndef VMS
  56. #include <fcntl.h>
  57. #endif    /* VMS */
  58. #endif /* OSK */
  59. #endif    /* AMIGA */
  60. #endif    /* DOS */
  61.  
  62. #if defined(ATARI) || defined(MTOS)
  63. #include <osbind.h>
  64. #include <stdio.h>
  65. #include <time.h>
  66. #define getchx() Crawcin()
  67. int kbhit(void);
  68. #endif
  69.  
  70. #define STANDARD    stderr    /* compatible with gnuplot philosophy */
  71.  
  72. #include "plot.h"
  73. #include "type.h"               /* own types */
  74. #include "matrix.h"
  75. #include "fit.h"
  76. #include "setshow.h"   /* for load_range */
  77.  
  78. enum marq_res { OK, ERROR, BETTER, WORSE };
  79. typedef enum marq_res marq_res_t;
  80.  
  81.  
  82. #define INFINITY    1e30
  83. #define NEARLY_ZERO 1e-30
  84. #define INITIAL_VALUE 1.0  /* create new variables with this value (was NEARLY_ZERO) */
  85. #define DELTA        0.001
  86. #define MAX_DATA    2048
  87. #define MAX_PARAMS  32
  88. #define MAX_LAMBDA  1e20
  89. #define MAX_VALUES_PER_LINE    32
  90. #define MAX_VARLEN  32
  91. #define START_LAMBDA 0.01
  92. #if defined(MSDOS) || defined(OS2) || defined(DOS386)
  93. #define PLUSMINUS   "\xF1"                      /* plusminus sign */
  94. #else
  95. #define PLUSMINUS   "+/-"
  96. #endif
  97. #define TEMPNAME "gfttemp.$$$"            /* used for msdos and os/2 */
  98.  
  99. static double epsilon       = 1e-5;        /* convergence criteria */
  100.  
  101. static char fit_script[127];
  102. static char logfile[128]    = "fit.log";
  103. static char *fit_index_def  = "FIT_INDEX";
  104. static char *FIXED        = "# FIXED";
  105. static char *GNUFITLOG        = "FIT_LOG";
  106. static char *FITLIMIT        = "FIT_LIMIT";
  107. static char *FITSCRIPT        = "FIT_SCRIPT";
  108. static char *DEFAULT_CMD    = "replot";         /* if no fitscript spec. */
  109.  
  110.  
  111. static FILE        *log_f = NULL;
  112.  
  113. static word        num_data,
  114.             num_params;
  115. static double        *fit_x,
  116.             *fit_y,
  117.             *err_data,
  118.             *a;
  119. static boolean        ctrlc_flag = FALSE;
  120.  
  121. static struct udft_entry func;
  122.  
  123. typedef char fixstr[MAX_VARLEN+1];
  124. static fixstr        *par_name;
  125.  
  126.  
  127. /*****************************************************************
  128.              internal Prototypes
  129. *****************************************************************/
  130.  
  131. static RETSIGTYPE ctrlc_handle __P((int an_int));
  132. static void ctrlc_setup __P((void));
  133. static marq_res_t marquardt __P((double a[], double **alpha, double *chisq,
  134.                     double *lambda, double *varianz));
  135. static boolean analyze __P((double a[], double **alpha, double beta[],
  136.                 double *chisq, double *varianz));
  137. static void calculate __P((double *yfunc, double **dyda, double a[]));
  138. static void call_gnuplot __P((double *par, double *data));
  139. static void show_fit __P((int i, double chisq, double last_chisq, double *a,
  140.                   double lambda, FILE *device));
  141. static boolean regress __P((double a[]));
  142. #if 0 /* not used */
  143. static int scan_num_entries __P((char *s, double vlist[]));
  144. #endif
  145. static boolean is_variable __P((char *s));
  146. static void copy_max __P((char *d, char *s, int n));
  147. static double getdvar __P((char *varname));
  148. static double createdvar __P((char *varname, double value));
  149. static void splitpath __P((char *s, char *p, char *f));
  150. #if defined(MSDOS) || defined(DOS386)
  151. static void subst __P((char *s, char from, char to));
  152. #endif
  153. static boolean fit_interrupt __P((void));
  154. static boolean is_empty __P((char *s));
  155.  
  156. /*****************************************************************
  157.     Useful macros
  158.     We avoid any use of varargs/stdargs (not good style but portable)
  159. *****************************************************************/
  160. #define Dblf(a)     {fprintf (STANDARD,a); fprintf (log_f,a);}
  161. #define Dblf2(a,b)    {fprintf (STANDARD,a,b); fprintf (log_f,a,b);}
  162. #define Dblf5(a,b,c,d,e) \
  163.         {fprintf (STANDARD,a,b,c,d,e); fprintf (log_f,a,b,c,d,e);}
  164.  
  165.  
  166. /*****************************************************************
  167.     This is called when a SIGINT occurs during fit
  168. *****************************************************************/
  169.  
  170. static RETSIGTYPE ctrlc_handle (an_int)
  171. int an_int;
  172. {
  173. #ifdef OS2
  174.     (void) signal(an_int, SIG_ACK);
  175. #else
  176.     /* reinstall signal handler (necessary on SysV) */
  177.     (void) signal(SIGINT, (sigfunc)ctrlc_handle);
  178. #endif
  179.     ctrlc_flag = TRUE;
  180. }
  181.  
  182.  
  183. /*****************************************************************
  184.     setup the ctrl_c signal handler
  185. *****************************************************************/
  186. static void ctrlc_setup()
  187. {
  188. /*
  189.  *  MSDOS defines signal(SIGINT) but doesn't handle it through
  190.  *  real interrupts. So there remain cases in which a ctrl-c may
  191.  *  be uncaught by signal. We must use kbhit() instead that really
  192.  *  serves the keyboard interrupt (or write an own interrupt func
  193.  *  which also generates #ifdefs)
  194.  *
  195.  *  I hope that other OSes do it better, if not... add #ifdefs :-(
  196.  */
  197. #if (defined(__EMX__) || !defined(MSDOS) && !defined(DOS386)) && !defined(ATARI) && !defined(MTOS)
  198.     (void) signal(SIGINT, (sigfunc)ctrlc_handle);
  199. #endif
  200. }
  201.  
  202.  
  203. /*****************************************************************
  204.     getch that handles also function keys etc.
  205. *****************************************************************/
  206. #if defined(MSDOS) || defined(DOS386)
  207. int getchx ()
  208. {
  209.     register int c = getch();
  210.     if ( !c  ||  c == 0xE0 ) {
  211.     c <<= 8;
  212.     c |= getch ();
  213.     }
  214.     return c;
  215. }
  216. #endif
  217.  
  218. /*****************************************************************
  219.     in case of fatal errors
  220. *****************************************************************/
  221. void error_ex ()
  222. {
  223.     char *sp;
  224.  
  225.     strncpy (fitbuf, "         ", 9);         /* start after GNUPLOT> */
  226.     sp = strchr (fitbuf, '\0');
  227.     while ( *--sp == '\n' )
  228.     ;
  229.     strcpy (sp+1, "\n\n");              /* terminate with exactly 2 newlines */
  230.     fprintf (STANDARD, fitbuf);
  231.       if ( log_f ) {
  232.     fprintf (log_f, "BREAK: %s", fitbuf);
  233.       fclose (log_f);
  234.       log_f = NULL;
  235.       }
  236.     if ( func.at ) {
  237.     free (func.at);             /* release perm. action table */
  238.     func.at = (struct at_type *) NULL;
  239.     }
  240.  
  241.     /* restore original SIGINT function */
  242. #if !defined(ATARI) && !defined(MTOS)
  243.     interrupt_setup ();
  244. #endif
  245.  
  246.     longjmp(env, TRUE);     /* bail out to command line */
  247. }
  248.  
  249.  
  250. /*****************************************************************
  251.     Marquardt's nonlinear least squares fit
  252. *****************************************************************/
  253. static marq_res_t marquardt (a, alpha, chisq, lambda, varianz)
  254. double a[];
  255. double **alpha;
  256. double *chisq;
  257. double *lambda;
  258. double *varianz;
  259. {
  260.     int         j;
  261.     static double   *da,        /* delta-step of the parameter */
  262.             *temp_a,        /* temptative new params set   */
  263.             **one_da,
  264.             *beta,
  265.             *tmp_beta,
  266.             **tmp_alpha,
  267.             **covar,
  268.             old_chisq;
  269.  
  270.     /* Initialization when lambda < 0 */
  271.  
  272.     if ( *lambda < 0 ) {    /* Get first chi-square check */
  273.     one_da        = matr (num_params, 1);
  274.     temp_a        = vec (num_params);
  275.     beta        = vec (num_params);
  276.     tmp_beta    = vec (num_params);
  277.     da        = vec (num_params);
  278.     covar        = matr (num_params, num_params);
  279.     tmp_alpha   = matr (num_params, num_params);
  280.  
  281.     *lambda  = -*lambda;      /* make lambda positive */
  282.     return analyze (a, alpha, beta, chisq, varianz) ? OK : ERROR;
  283.     }
  284.     old_chisq = *chisq;
  285.  
  286.     for ( j=0 ; j<num_params ; j++ ) {
  287.     memcpy (covar[j], alpha[j], num_params*sizeof(double));
  288.     covar[j][j] *= 1 + *lambda;
  289.     one_da[j][0] = beta [j];
  290.     }
  291.  
  292.     solve (covar, num_params, one_da, 1);    /* Equation solution */
  293.  
  294.     for ( j=0 ; j<num_params ; j++ )
  295.     da[j] = one_da[j][0];            /* changes in paramss */
  296.  
  297.     /* once converged, free dynamic allocated vars */
  298.  
  299.     if ( *lambda == 0.0 ) {
  300.     free (beta);
  301.     free (tmp_beta);
  302.     free (da);
  303.     free (temp_a);
  304.     free_matr (one_da, num_params);
  305.     free_matr (tmp_alpha, num_params);
  306.     free_matr (covar, num_params);
  307.         return OK;
  308.     }
  309.  
  310.     /* check if trial did ameliorate sum of squares */
  311.  
  312.     for ( j=0 ; j<num_params ; j++ ) {
  313.     temp_a[j] = a[j] + da[j];
  314.     tmp_beta[j] = beta [j];
  315.     memcpy (tmp_alpha[j], alpha[j], num_params*sizeof(double));
  316.     }
  317.  
  318.     if ( !analyze (temp_a, tmp_alpha, tmp_beta, chisq, varianz) )
  319.     return ERROR;
  320.  
  321.     if ( *chisq < old_chisq ) {     /* Success, accept new solution */
  322.     if ( *lambda > 1e-20 )
  323.         *lambda /= 10;
  324.     old_chisq = *chisq;
  325.     for ( j=0 ; j<num_params ; j++ ) {
  326.         memcpy (alpha[j], tmp_alpha[j], num_params*sizeof(double));
  327.         beta[j] = tmp_beta[j];
  328.         a[j] = temp_a[j];
  329.     }
  330.     return BETTER;
  331.     }
  332.     else {                /* failure, increase lambda and return */
  333.     *lambda *= 10;
  334.     *chisq = old_chisq;
  335.     return WORSE;
  336.     }
  337. }
  338.  
  339.  
  340. /*****************************************************************
  341.     compute chi-square and numeric derivations
  342. *****************************************************************/
  343. static boolean analyze (a, alpha, beta, chisq, varianz)
  344. double a[];
  345. double **alpha;
  346. double beta[];
  347. double *chisq;
  348. double *varianz;
  349. {
  350.  
  351. /*
  352.  *  used by marquardt to evaluate the linearized fitting matrix alpha
  353.  *  and vector beta
  354.  */
  355.  
  356.     word    k, j, i;
  357.     double  wt, sig2i, dy, **dyda, *yfunc;
  358. #if !defined(ATARI) && !defined(MTOS)
  359.     double *edp, *yp, *fp, **dr, **ar, *dc, *bp, *ac;
  360. #endif
  361.  
  362.     yfunc = vec (num_data);
  363.     dyda = matr (num_data, num_params);
  364.  
  365.     /* initialize alpha beta */
  366.     for ( j=0 ; j<num_params ; j++ ) {
  367.     for ( k=0 ; k<=j ; k++ )
  368.         alpha[j][k] = 0;
  369.     beta[j] = 0;
  370.     }
  371.  
  372.     *chisq = *varianz = 0;
  373.     calculate (yfunc, dyda, a);
  374.  
  375. #if !defined(ATARI) && !defined(MTOS)
  376.     edp = err_data;
  377.     yp = fit_y;
  378.     fp = yfunc;
  379.     dr = dyda;
  380.  
  381.     /* Summation loop over all data */
  382.     for ( i=0 ; i<num_data ; i++, dr++ ) {
  383.     /* The original code read: sig2i = 1/(*edp * *edp++); */
  384.     /* There are some compilers that evaluate the operation */
  385.     /* from right to left, although it is an error to do so. */
  386.     /* Hence the following modification: */
  387.     sig2i = 1/(*edp * *edp);
  388.     edp++;
  389.     dy = *yp++ - *fp++;
  390.     *varianz += dy*dy;
  391.     ar = alpha;
  392.         dc = *dr;
  393.     bp = beta;
  394.     for ( j=0 ; j<num_params ; j++ ) {
  395.         wt = *dc++ * sig2i;
  396.         ac = *ar++;
  397.             for ( k=0 ; k<=j ; k++ )
  398.         *ac++ += wt * (*dr)[k];
  399.         *bp++ += dy * wt;
  400.     }
  401.     *chisq += dy*dy*sig2i;            /* and find chi-square */
  402.     }
  403. #else
  404.      /* Summation loop over all data */
  405.     for ( i=0 ; i<num_data ; i++) {
  406.         sig2i = 1/(err_data[i] * err_data[i]);
  407.         dy = fit_y[i] - yfunc[i];
  408.     *varianz += dy*dy;
  409.         for ( j=0 ; j<num_params ; j++ ) {
  410.         wt = dyda[i][j] * sig2i;
  411.             for ( k=0 ; k<=j ; k++ )
  412.         alpha[j][k] += wt * dyda[i][k];
  413.             beta[j] += dy * wt;
  414.     }
  415.     *chisq += dy*dy*sig2i;            /* and find chi-square */
  416.     }
  417. #endif
  418.     *varianz /= num_data;
  419.     for ( j=1 ; j<num_params ; j++ )        /* fill in the symmetric side */
  420.     for ( k=0 ; k<=j-1 ; k++ )
  421.         alpha[k][j] = alpha [j][k];
  422.     free (yfunc);
  423.     free_matr (dyda, num_data);
  424.     return TRUE;
  425. }
  426.  
  427.  
  428. /*****************************************************************
  429.     compute function values and partial derivatives of chi-square
  430. *****************************************************************/
  431. static void calculate (yfunc, dyda, a)
  432. double *yfunc;
  433. double **dyda;
  434. double a[];
  435. {
  436.     word    k, p;
  437.     double  tmp_a;
  438.     double  *tmp_low,
  439.         *tmp_high,
  440.         *tmp_pars;
  441.  
  442.     tmp_low = vec (num_data);          /* numeric derivations */
  443.     tmp_high = vec (num_data);
  444.     tmp_pars = vec (num_params);
  445.  
  446.     /* first function values */
  447.  
  448.     call_gnuplot (a, yfunc);
  449.  
  450.     /* then derivatives */
  451.  
  452.     for ( p=0 ; p<num_params ; p++ )
  453.     tmp_pars[p] = a[p];
  454.     for ( p=0 ; p<num_params ; p++ ) {
  455.     tmp_pars[p] = tmp_a = fabs(a[p]) < NEARLY_ZERO ? NEARLY_ZERO : a[p];
  456. /*
  457.  *  the more exact method costs double execution time and is therefore
  458.  *  commented out here. Change if you like!
  459.  */
  460. /*      tmp_pars[p] *= 1-DELTA;
  461.     call_gnuplot (tmp_pars, tmp_low);  */
  462.  
  463.     tmp_pars[p] = tmp_a * (1+DELTA);
  464.     call_gnuplot (tmp_pars, tmp_high);
  465.     for ( k=0 ; k<num_data ; k++ )
  466.  
  467. /*            dyda[k][p] = (tmp_high[k] - tmp_low[k])/(2*tmp_a*DELTA); */
  468.  
  469.         dyda[k][p] = (tmp_high[k] - yfunc[k])/(tmp_a*DELTA);
  470.         tmp_pars[p] = a[p];
  471.     }
  472.  
  473.     free (tmp_low);
  474.     free (tmp_high);
  475.     free (tmp_pars);
  476. }
  477.  
  478.  
  479. /*****************************************************************
  480.     call internal gnuplot functions
  481. *****************************************************************/
  482. static void call_gnuplot (par, data)
  483. double *par;
  484. double *data;
  485. {
  486.     word    i;
  487.     struct value v;
  488.  
  489.     /* set parameters first */
  490.     for ( i=0 ; i<num_params ; i++ ) {
  491.     Gcomplex (&v, par[i], 0.0);
  492.     setvar (par_name[i], v);
  493.     }
  494.  
  495.     for (i = 0; i < num_data; i++) {
  496.     /* fit_index for multi-branch functions */
  497.     (void) Ginteger (&func.dummy_values[0], (i+1));
  498.     setvar (fit_index, func.dummy_values[0]);
  499.     /* calculate fit-function value */
  500.     (void) Gcomplex (&func.dummy_values[0], fit_x[i], 0.0);
  501.     evaluate_at (func.at,&v);
  502.     if ( undefined )
  503.         Eex ("Undefined value during function evaluation")
  504.     data[i] = real(&v);
  505.     }
  506.     (void) Ginteger (&func.dummy_values[0], 0);
  507.     setvar (fit_index, func.dummy_values[0]);
  508. }
  509.  
  510.  
  511. /*****************************************************************
  512.     handle user interrupts during fit
  513. *****************************************************************/
  514. static boolean fit_interrupt ()
  515. {
  516.     char c;
  517.  
  518.     while ( TRUE ) {
  519.     fprintf (STANDARD, "\n\n(S)top fit, (C)ontinue, (E)xecute:  ");
  520.         /* toupper is a macro that might call fgetc twice on the Amiga */
  521.         c = fgetc (stdin);
  522.     switch (toupper(c)) {
  523.     case 'S' :  fprintf (STANDARD, "Stop.");
  524.             return FALSE;
  525.     case 'C' :  fprintf (STANDARD, "Continue.");
  526.             return TRUE;
  527.     case 'E' :  {
  528.             int i;
  529.             struct value v;
  530.             char *tmp;
  531.  
  532.             tmp = *fit_script ? fit_script : DEFAULT_CMD;
  533.             fprintf (STANDARD, "executing: %s", tmp);
  534.             /* set parameters visible to gnuplot */
  535.             for ( i=0 ; i<num_params ; i++ ) {
  536.                 Gcomplex (&v, a[i], 0.0);
  537.                 setvar (par_name[i], v);
  538.             }
  539.             sprintf (input_line, tmp);
  540.             (void) do_line ();
  541.             }
  542.     }
  543.     }
  544. }
  545.  
  546.  
  547. /*****************************************************************
  548.     frame routine for the marquardt-fit
  549. *****************************************************************/
  550. static boolean regress (a)
  551. double a[];
  552. {
  553.     double    **covar,
  554.         **correl,
  555.         *dpar,
  556.         **alpha,
  557.         chisq,
  558.         last_chisq,
  559.         varianz,
  560.         lambda;
  561.     word    i, j;
  562.     marq_res_t    res;
  563.     struct value val;
  564.  
  565.     chisq   = last_chisq = INFINITY;
  566.     alpha   = matr (num_params, num_params);
  567.     lambda  = -START_LAMBDA;            /* use sign as flag */
  568.     i = 0;                    /* iteration counter  */
  569.  
  570.     /* ctrlc now serves as Hotkey */
  571.     ctrlc_setup ();
  572.  
  573.     /* Initialize internal variables and 1st chi-square check */
  574.     if ( (res = marquardt (a, alpha, &chisq, &lambda, &varianz)) == ERROR )
  575.     Eex ("FIT: error occured during fit")
  576.     res = BETTER;
  577.  
  578.     Dblf ("\nInitial set of free parameters:\n")
  579.     show_fit (i, chisq, chisq, a, lambda, STANDARD);
  580.     show_fit (i, chisq, chisq, a, lambda, log_f);
  581.  
  582.     /* MAIN FIT LOOP: do the regression iteration */
  583.  
  584.     do {
  585. /*
  586.  *  MSDOS defines signal(SIGINT) but doesn't handle it through
  587.  *  real interrupts. So there remain cases in which a ctrl-c may
  588.  *  be uncaught by signal. We must use kbhit() instead that really
  589.  *  serves the keyboard interrupt (or write an own interrupt func
  590.  *  which also generates #ifdefs)
  591.  *
  592.  *  I hope that other OSes do it better, if not... add #ifdefs :-(
  593.  *  EMX does not have kbhit.
  594.  */
  595. #if ((defined(MSDOS) || defined(DOS386)) && !defined(__EMX__)) || defined(ATARI) || defined(MTOS)
  596.      if ( kbhit () ) {
  597.         do {getchx ();} while ( kbhit() );
  598.         ctrlc_flag = TRUE;
  599.     }
  600. #endif
  601.  
  602.         if ( ctrlc_flag ) {
  603.         show_fit (i, chisq, last_chisq, a, lambda, STANDARD);
  604.         ctrlc_flag = FALSE;
  605.         if ( !fit_interrupt () )           /* handle keys */
  606.         break;
  607.     }
  608.     if ( res == BETTER ) {
  609.             i++;
  610.             last_chisq = chisq;
  611.     }
  612.         if ( (res = marquardt (a, alpha, &chisq, &lambda, &varianz)) == BETTER )
  613.         show_fit (i, chisq, last_chisq, a, lambda, STANDARD);
  614.     }
  615.     while ( res != ERROR  &&  lambda < MAX_LAMBDA  &&
  616.         (res == WORSE  ||  ( chisq > NEARLY_ZERO ? (last_chisq-chisq)/chisq : (last_chisq-chisq)) > epsilon) );
  617.  
  618.     /* restore original SIGINT function */
  619. #if !defined(ATARI) && !defined(MTOS)
  620.     interrupt_setup ();
  621. #endif
  622.     Dblf2 ("\nAfter %d iterations the fit converged.\n", i)
  623.     Dblf2 ("final sum of squares residuals    : %g\n", chisq)
  624.     if (chisq > NEARLY_ZERO) {
  625.         Dblf2 ("rel. change during last iteration : %g\n\n", (chisq-last_chisq)/chisq)
  626.      } else {
  627.         Dblf2 ("abs. change during last iteration : %g\n\n", (chisq-last_chisq))
  628.      }
  629.  
  630.     if ( res == ERROR )
  631.     Eex ("FIT: error occured during fit")
  632.  
  633.     /* compute errors in the parameters */
  634.  
  635.     if (varianz > NEARLY_ZERO) { /* what should we do if varianz=0 (exact fit) ? */
  636.     for ( i=0 ; i<num_params; i++ )
  637.     for ( j=0 ; j<num_params; j++ )
  638.         alpha[i][j] /= varianz;
  639.      }
  640.      
  641.     /* get covariance-, Korrelations- and Kurvature-Matrix */
  642.     /* and errors in the parameters              */
  643.  
  644.     covar   = matr (num_params, num_params);
  645.     correl  = matr (num_params, num_params);
  646.     dpar = vec (num_params);
  647.  
  648.     inverse (alpha, covar, num_params);
  649.  
  650.     Dblf ("Final set of parameters            68.3%% confidence interval (one at a time)\n")
  651.     Dblf ("=======================            =========================================\n\n")
  652.     for ( i=0 ; i<num_params; i++ ) {
  653.     dpar[i] = sqrt (covar[i][i]);
  654.     Dblf5 ("%-15.15s = %-15g  %-3.3s %g\n", par_name[i], a[i], PLUSMINUS, dpar[i])
  655.     }
  656.  
  657.     Dblf ("\n\ncorrelation matrix of the fit parameters:\n\n")
  658.     Dblf ("               ")
  659.     for ( j=0 ; j<num_params ; j++ )
  660.     Dblf2 ("%-6.6s ", par_name[j])
  661.     Dblf ("\n")
  662.  
  663.     for ( i=0 ; i<num_params; i++ ) {
  664.     Dblf2 ("%-15.15s", par_name[i])
  665.     for ( j=0 ; j<num_params; j++ ) {
  666.         correl[i][j] = covar[i][j] / (dpar[i]*dpar[j]);
  667.         Dblf2 ("%6.3f ", correl[i][j])
  668.         }
  669.     Dblf ("\n")
  670.     }
  671.  
  672.     /* restore last parameter's value (not done by calculate) */
  673.     Gcomplex (&val, a[num_params-1], 0.0);
  674.     setvar (par_name[num_params-1], val);
  675.  
  676.     /* call destruktor for allocated vars */
  677.     lambda = 0;
  678.     marquardt (a, alpha, &chisq, &lambda, &varianz);
  679.  
  680.     free_matr (covar, num_params);
  681.     free_matr (alpha, num_params);
  682.     free_matr (correl, num_params);
  683.     free (dpar);
  684.  
  685.     return TRUE;
  686. }
  687.  
  688.  
  689. /*****************************************************************
  690.     display actual state of the fit
  691. *****************************************************************/
  692. static void show_fit (i, chisq, last_chisq, a, lambda, device)
  693. int i;
  694. double chisq;
  695. double last_chisq;
  696. double *a;
  697. double lambda;
  698. FILE *device;
  699. {
  700.     int k;
  701.  
  702.     fprintf (device, "\n\n\
  703. Iteration %d\n\
  704. chisquare : %-15g   relative deltachi2 : %g\n\
  705. deltachi2 : %-15g   limit for stopping : %g\n\
  706. lambda      : %g\n\nactual set of parameters\n\n",
  707. i, chisq, chisq > NEARLY_ZERO ? (chisq - last_chisq)/chisq : 0.0, chisq - last_chisq, epsilon, lambda);
  708.     for ( k=0 ; k<num_params ; k++ )
  709.     fprintf (device, "%-15.15s = %g\n", par_name[k], a[k]);
  710. }
  711.  
  712.  
  713.  
  714. /*****************************************************************
  715.     is_empty: check for valid string entries
  716. *****************************************************************/
  717. static boolean is_empty (s)
  718. char *s;
  719. {
  720.     while ( *s == ' '  ||  *s == '\t'  ||  *s == '\n' )
  721.     s++;
  722.     return (boolean)( *s == '#'  ||  *s == '\0' );
  723. }
  724.  
  725.  
  726. /*****************************************************************
  727.     get next word of a multi-word string, advance pointer
  728. *****************************************************************/
  729. char *get_next_word (s, subst)
  730. char **s;
  731. char *subst;
  732. {
  733.     char *tmp = *s;
  734.  
  735.     while ( *tmp==' ' || *tmp=='\t' || *tmp=='=' )
  736.     tmp++;
  737.     if ( *tmp=='\n' || *tmp=='\0' )                    /* not found */
  738.     return NULL;
  739.     if ( (*s = strpbrk (tmp, " =\t\n")) == NULL )
  740.     *s = tmp + strlen(tmp);
  741.     *subst = **s;
  742.     *(*s)++ = '\0';
  743.     return tmp;
  744. }
  745.  
  746.  
  747.  
  748. /*****************************************************************
  749.     get valid numerical entries
  750. *****************************************************************/
  751. #if 0 /* not used */
  752. static int scan_num_entries (s, vlist)
  753. char *s;
  754. double vlist[];
  755. {
  756.     int num = 0;
  757.     char c, *tmp;
  758.  
  759.     while ( (tmp = get_next_word (&s, &c)) != NULL
  760.                     &&  num <= MAX_VALUES_PER_LINE )
  761.     if ( !sscanf (tmp, "%lf", &vlist[++num]) )
  762.         Eex ("syntax error in data file")
  763.     return num;
  764. }
  765. #endif
  766.  
  767. /*****************************************************************
  768.     check for variable identifiers
  769. *****************************************************************/
  770. static boolean is_variable (s)
  771. char *s;
  772. {
  773.     while ( *s ) {
  774.     if ( !isalnum(*s) && *s!='_' )
  775.         return FALSE;
  776.     s++;
  777.     }
  778.     return TRUE;
  779. }
  780.  
  781.  
  782. /*****************************************************************
  783.     strcpy but max n chars
  784. *****************************************************************/
  785. static void copy_max (d, s, n)
  786. char *d;
  787. char *s;
  788. int n;
  789. {
  790.     strncpy (d, s, n);
  791.     if ( strlen(s) >= n )
  792.     d[n-1] = '\0';
  793. }
  794.  
  795.  
  796. /*****************************************************************
  797.     portable implementation of strnicmp (hopefully)
  798. *****************************************************************/
  799.  
  800. #ifdef NEED_STRNICMP
  801. int strnicmp (s1, s2, n)
  802. char *s1;
  803. char *s2;
  804. int n;
  805. {
  806.     char c1,c2;
  807.  
  808.     if(n==0) return 0;
  809.  
  810.     do {
  811.     c1 = *s1++; if(islower(c1)) c1=toupper(c1);
  812.     c2 = *s2++; if(islower(c2)) c2=toupper(c2);
  813.     } while(c1==c2 && c1 && c2 && --n>0);
  814.  
  815.     if(n==0 || c1==c2) return 0;
  816.     if(c1=='\0' || c1<c2) return 1;
  817.     return -1;
  818. }
  819. #endif
  820.  
  821.  
  822. /*****************************************************************
  823.     first time settings
  824. *****************************************************************/
  825. void init_fit ()
  826. {
  827.     struct value val;
  828.  
  829.     /* index used for multi-branch functions */
  830.     fit_index = fit_index_def;
  831.     val.type = INTGR;
  832.     val.v.int_val = 0;
  833.     setvar (fit_index, val);
  834.     func.at = (struct at_type *) NULL;      /* need to parse 1 time */
  835. }
  836.  
  837.  
  838. /*****************************************************************
  839.         Set a GNUPLOT user-defined variable
  840. ******************************************************************/
  841.  
  842. void setvar (varname, data)
  843. char *varname;
  844. struct value data;
  845. {
  846.     register struct udvt_entry *udv_ptr = first_udv,
  847.                    *last = first_udv;
  848.  
  849.     /* check if it's already in the table... */
  850.  
  851.     while (udv_ptr) {
  852.         last = udv_ptr;
  853.         if ( !strcmp (varname, udv_ptr->udv_name))
  854.         break;
  855.         udv_ptr = udv_ptr->next_udv;
  856.     }
  857.  
  858.     if ( !udv_ptr ) {                /* generate new entry */
  859.     last->next_udv = udv_ptr = (struct udvt_entry *)
  860.       alloc ((unsigned int)sizeof(struct udvt_entry), "setvar");
  861.     udv_ptr->next_udv = NULL;
  862.     }
  863.     copy_max (udv_ptr->udv_name, varname, sizeof(udv_ptr->udv_name));
  864.     udv_ptr->udv_value = data;
  865.     udv_ptr->udv_undef = FALSE;
  866. }
  867.  
  868.  
  869.  
  870. /*****************************************************************
  871.     Read INTGR Variable value, return 0 if undefined or wrong type
  872. *****************************************************************/
  873. int getivar (varname)
  874. char *varname;
  875. {
  876.     register struct udvt_entry *udv_ptr = first_udv;
  877.  
  878.     while (udv_ptr) {
  879.         if ( !strcmp (varname, udv_ptr->udv_name))
  880.         return udv_ptr->udv_value.type == INTGR
  881.                ? udv_ptr->udv_value.v.int_val    /* valid */
  882.                : 0;                /* wrong type */
  883.             udv_ptr = udv_ptr->next_udv;
  884.     }
  885.     return 0;                        /* not in table */
  886. }
  887.  
  888.  
  889. /*****************************************************************
  890.     Read DOUBLE Variable value, return 0 if undefined or wrong type
  891.    I dont think it's a problem that it's an integer - div
  892. *****************************************************************/
  893. static double getdvar (varname)
  894. char *varname;
  895. {
  896.     register struct udvt_entry *udv_ptr = first_udv;
  897.  
  898.     for ( ; udv_ptr ; udv_ptr=udv_ptr->next_udv)
  899.         if ( strcmp (varname, udv_ptr->udv_name)==0)
  900.           return real(&(udv_ptr->udv_value));
  901.  
  902.     /* get here => not found */
  903.     return 0;
  904. }
  905.  
  906. /* like getdvar, but
  907.  * - convert it from integer to real if necessary
  908.  * - create it with value INITIAL_VALUE if not found or undefined
  909.  */
  910. static double createdvar(varname,value)
  911. char *varname;
  912. double value;
  913. {
  914.     register struct udvt_entry *udv_ptr = first_udv;
  915.  
  916.     for ( ; udv_ptr ; udv_ptr=udv_ptr->next_udv)
  917.         if ( strcmp (varname, udv_ptr->udv_name)==0) {
  918.             if (udv_ptr->udv_undef) {
  919.                 udv_ptr->udv_undef=0;
  920.                 Gcomplex(&udv_ptr->udv_value, value, 0.0);
  921.             } else if (udv_ptr->udv_value.type==INTGR) {
  922.                 Gcomplex(&udv_ptr->udv_value, (double)udv_ptr->udv_value.v.int_val, 0.0);
  923.             }
  924.            return real(&(udv_ptr->udv_value));
  925.          }
  926.  
  927.     /* get here => not found */
  928.  
  929.     {
  930.         struct value a;
  931.         Gcomplex(&a, value, 0.0);
  932.         setvar(varname, a);
  933.     }
  934.     
  935.     return value;
  936. }
  937.     
  938.  
  939. /*****************************************************************
  940.     Split Identifier into path and filename
  941. *****************************************************************/
  942. static void splitpath (s, p, f)
  943. char *s;
  944. char *p;
  945. char *f;
  946. {
  947.     register char *tmp;
  948.     tmp = s + strlen(s) - 1;
  949.     while ( *tmp != '\\'  &&  *tmp != '/'  &&  *tmp != ':'  &&  tmp-s >= 0 )
  950.     tmp--;
  951.     strcpy (f, tmp+1);
  952.     strncpy (p, s, tmp-s+1);
  953.     p[tmp-s+1] = '\0';
  954. }
  955.  
  956.  
  957. #if defined(MSDOS) || defined(DOS386)
  958. /*****************************************************************
  959.     Character substitution
  960. *****************************************************************/
  961. #ifdef ANSI_C
  962. static void subst (char *s, char from, char to)
  963. #else
  964. static void subst (s, from, to)
  965. char *s;
  966. char from;
  967. char to;
  968. #endif
  969. {
  970.     while ( *s ) {
  971.     if ( *s == from )
  972.         *s = to;
  973.     s++;
  974.     }
  975. }
  976. #endif
  977.  
  978.  
  979. /*****************************************************************
  980.     write the actual parameters to start parameter file
  981. *****************************************************************/
  982. void update (pfile)
  983. char *pfile;
  984. {
  985.     char    fnam[256],
  986.         path[256],
  987.         tmpname[256],
  988.         sstr[256],
  989.         pname[64],
  990.         tail[127],
  991.         *s = sstr,
  992. #if !defined(AMIGA_SC_6_1)
  993.         *tmpn,
  994. #endif
  995.         *tmp,
  996.         c;
  997.     FILE        *of,
  998.         *nf;
  999.     double    pval;
  1000.  
  1001.     /* split into path and filename */
  1002.  
  1003. #if defined(MSDOS) || defined(DOS386)
  1004.     /* don't get problems during system() call */
  1005.     subst (pfile, '/', '\\');
  1006. #endif
  1007.     splitpath (pfile, path, fnam);
  1008.     if ( !(of = fopen (pfile, "r")) )
  1009.     Eex2 ("parameter file %s could not be read", pfile)
  1010. #if defined(AMIGA_SC_6_1)
  1011.     sprintf (tmpname, "%sgf%08lx", path, FindTask (NULL));
  1012. #else /* !AMIGA_SC_6_1 */
  1013.     /* Under MSDOS rename will not work properly if TMP points to another
  1014.      * drive/dir. In this case use 'path' directory for tempfile
  1015.      */
  1016. #if defined (MSDOS)  ||  defined (OS2)
  1017.     sprintf (tmpname, "%s%s", path, TEMPNAME);
  1018. #else
  1019. #ifdef HAVE_TEMPNAM
  1020.     tmpn = tempnam (path, "gnuft");
  1021. #else
  1022.     tmpn = tmpnam(NULL);
  1023. #endif
  1024.  
  1025.     if ( tmpn == NULL )
  1026.     Eex ("no more unique names possible during updating")
  1027.     strcpy (tmpname, tmpn);
  1028.     free (tmpn);
  1029. #endif
  1030. #endif /* !AMIGA_SC_6_1 */
  1031.     if ( !(nf = fopen (tmpname, "w")) )
  1032.     Eex ("could not open temporary file")
  1033.     while ( TRUE ) {
  1034.     if ( !fgets (s = sstr, sizeof(sstr), of) )    /* EOF found */
  1035.         break;
  1036.     if ( is_empty(s) ) {
  1037.         fprintf (nf, s);                /* preserve comments */
  1038.             continue;
  1039.         }
  1040.         if ( (tmp = strchr (s, '#')) != NULL ) {
  1041.         copy_max (tail, tmp, sizeof(tail));
  1042.         *tmp = '\0';
  1043.     }
  1044.     else
  1045.         strcpy (tail, "\n");
  1046.     tmp = get_next_word (&s, &c);
  1047.     if ( !is_variable (tmp)  ||  strlen(tmp) > MAX_VARLEN ) {
  1048.         fclose (of);
  1049.         Eex2 ("syntax error in parameter file %s", fnam)
  1050.     }
  1051.     copy_max (pname, tmp, sizeof(pname));
  1052.     /* next must be '=' */
  1053.     if ( c != '=' ) {
  1054.         tmp = strchr (s, '=');
  1055.         if ( tmp == NULL ) {
  1056.         fclose (of);
  1057.         Eex2 ("syntax error in parameter file %s", fnam)
  1058.         }
  1059.         s = tmp+1;
  1060.     }
  1061.     tmp = get_next_word (&s, &c);
  1062.     if ( !sscanf (tmp, "%lf", &pval) ) {
  1063.         fclose (of);
  1064.         Eex2 ("syntax error in parameter file %s", fnam)
  1065.     }
  1066.     if ( (tmp = get_next_word (&s, &c)) != NULL ) {
  1067.         fclose (of);
  1068.         Eex2 ("syntax error in parameter file %s", fnam)
  1069.     }
  1070.  
  1071.         /* now modify */
  1072.  
  1073.     if ( (pval = getdvar (pname)) == 0 )
  1074.         pval = (double) getivar (pname);
  1075.  
  1076.     sprintf (sstr, "%g", pval);
  1077.     if ( !strchr (sstr, '.')  &&  !strchr (sstr, 'e') )
  1078.         strcat (sstr, ".0");                /* assure CMPLX-type */
  1079.     fprintf (nf, "%-15.15s = %-15.15s   %s", pname, sstr, tail);
  1080.     }
  1081.     if ( fclose (of)  ||  fclose (nf) )
  1082.     Eex ("I/O error during update")
  1083.  
  1084. #ifdef VMS
  1085.     if ( delete (pfile) )     /* implies write permission !!!! */
  1086. #else    /* VMS */
  1087.     if ( unlink (pfile) )     /* implies write permission !!!! */
  1088. #endif    /* VMS */
  1089.       Eex2 ("parameter file %s could not be deleted", pfile)
  1090.  
  1091. #if defined(MSDOS) || defined(OS2)
  1092.     sprintf (sstr, "ren %s %s", tmpname, fnam);
  1093. #ifdef _Windows
  1094.     if ( winsystem (sstr) )
  1095. #else
  1096.     if ( system (sstr) )
  1097. #endif
  1098. #else                   /* implies write permission !!!! */
  1099.     if ( rename (tmpname, pfile) )
  1100. #endif
  1101.       Eex3 ("unable to rename %s to %s", tmpname, pfile)
  1102. }
  1103.  
  1104.  
  1105.  
  1106. /*****************************************************************
  1107.     Interface to the classic gnuplot-software
  1108. *****************************************************************/
  1109.  
  1110. void do_fit ()
  1111. {
  1112.     int autorange=3;  /* yes */
  1113.     double min, max;  /* range to fit */
  1114.     int dummy_token=-1; /* local variable name */
  1115.  
  1116.     int     i;
  1117.     int columns;
  1118.     double    v[3];
  1119.     double tmpd;
  1120.     time_t    timer;
  1121.     int token1, token2, token3;
  1122.     char *tmp;
  1123.     extern struct udft_entry *dummy_func;
  1124.     extern char dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
  1125.     extern char c_dummy_var[MAX_NUM_VAR][MAX_ID_LEN+1];
  1126.     extern int c_token;
  1127.    extern int df_line_number;
  1128.    
  1129. /* first look for a restricted fit range */
  1130.  
  1131.     if (equals(c_token, "[")) {
  1132.     c_token++;
  1133.     if (isletter(c_token)) {
  1134.         if (equals(c_token + 1, "=")) {
  1135.         dummy_token = c_token;
  1136.         c_token += 2;
  1137.         }
  1138.             /* else parse it as a xmin expression */
  1139.     }
  1140.     autorange = load_range(FIRST_X_AXIS,&min, &max, autorange);
  1141.     if (!equals(c_token, "]"))
  1142.         int_error("']' expected", c_token);
  1143.     c_token++;
  1144.     }
  1145.  
  1146.  
  1147. /* now compile the function */
  1148.  
  1149.     token1 = c_token;
  1150.  
  1151.     if (func.at) {
  1152.     free(func.at);
  1153.     func.at=NULL;  /* in case perm_at() does int_error */
  1154.     }
  1155.  
  1156.     dummy_func = &func;
  1157.  
  1158.     if (dummy_token >= 0)
  1159.     copy_str(c_dummy_var[0], dummy_token, MAX_ID_LEN);
  1160.     else
  1161.     strcpy(c_dummy_var[0], dummy_var[0]);
  1162.     
  1163.     func.at = perm_at();
  1164.     dummy_func = NULL;
  1165.  
  1166.     token2 = c_token;
  1167.  
  1168. /* use datafile module to parse the datafile and qualifiers */
  1169.  
  1170.     columns = df_open(3);  /* up to 3 using specs allowed */
  1171.     if (columns==1)
  1172.     int_error("Need 2 or 3 using specs", c_token);
  1173.     /* defer actually reading the data until we have parsed the rest of the line */
  1174.     
  1175.     token3 = c_token;
  1176.  
  1177.     tmpd = getdvar (FITLIMIT);          /* get epsilon if given explicitly */
  1178.     if ( tmpd < 1.0  &&  tmpd > 0.0 )
  1179.         epsilon = tmpd;
  1180.  
  1181.     *fit_script = '\0';                   /* FIT_SKIP 1 means each 2nd value */
  1182.     if ( (tmp = getenv (FITSCRIPT)) != NULL )
  1183.     strcpy (fit_script, tmp);
  1184.  
  1185.     tmp = getenv (GNUFITLOG);          /* open logfile */
  1186.     if ( tmp != NULL ) {
  1187.         char *tmp2 = &tmp[strlen(tmp)-1];
  1188.         if ( *tmp2 == '/'  ||  *tmp2 == '\\' )
  1189.             sprintf (logfile, "%s%s", tmp, logfile);
  1190.         else
  1191.             strcpy (logfile, tmp);
  1192.     }
  1193.     if ( !log_f && /* div */ !(log_f = fopen (logfile, "a")) )
  1194.     Eex2 ("could not open log-file %s", logfile)
  1195.     fprintf (log_f, "\n\n*******************************************************************************\n");
  1196.     time (&timer);
  1197.     fprintf (log_f, "%s\n\n", ctime (&timer));
  1198.     {    char line[MAX_LINE_LEN];
  1199.         capture(line, token2,token3-1,MAX_LINE_LEN);
  1200.             fprintf (log_f, "FIT:    data read from %s\n", line);
  1201.     }
  1202.  
  1203.     fit_x = vec (MAX_DATA);         /* start with max. value */
  1204.     fit_y = vec (MAX_DATA);
  1205.  
  1206. /* first read in experimental data */
  1207.  
  1208.     err_data = vec (MAX_DATA);
  1209.     num_data = 0;
  1210.  
  1211.     while ( (i=df_readline(v, 3)) != EOF ) {
  1212.  
  1213.         if (i==0)
  1214.             continue;
  1215.  
  1216.         if ( !(autorange&1) && (v[0] < min))
  1217.             continue;  /* skip this one - div */
  1218.         if ( !(autorange&2) && (v[0] > max))
  1219.             continue;  /* and again */
  1220.  
  1221.         if ( num_data==MAX_DATA ) {
  1222.             df_close();
  1223.             Eex2 ("max. # of datapoints %d exceeded", MAX_DATA);
  1224.         }
  1225.  
  1226.         switch (i)
  1227.         {
  1228.             case DF_UNDEFINED:
  1229.             case DF_FIRST_BLANK:
  1230.             case DF_SECOND_BLANK:
  1231.                 continue;
  1232.             case 0:
  1233.                 Eex2("bad data on line %d", df_line_number);
  1234.                 break;
  1235.             case 1:
  1236.                 v[1]=v[0];
  1237.                 v[0]=df_line_number;
  1238.                 break;
  1239.         }
  1240.                 
  1241.         fit_x[num_data] = v[0];
  1242.         fit_y[num_data] = v[1];
  1243.         /* datafile module is such that if columns==3, it
  1244.          * will only return lines with 3 data read
  1245.          */
  1246.         err_data[num_data++] = (columns==3 ? v[2] : 1);
  1247.     }
  1248.     df_close();
  1249.  
  1250.     /* now resize fields */
  1251.     redim_vec (&fit_x, num_data);
  1252.     redim_vec (&fit_y, num_data);
  1253.     redim_vec (&err_data, num_data);
  1254.     
  1255.     fprintf (log_f, "        #datapoints = %d\n", num_data);
  1256.     if ( columns==2 )
  1257.     fprintf (log_f, "        y-errors assumed equally distributed\n\n");
  1258.  
  1259.     {    char line[MAX_LINE_LEN];
  1260.         capture(line, token1, token2-1,MAX_LINE_LEN);        
  1261.         fprintf (log_f, "function used for fitting: %s\n", line);
  1262.     }
  1263.  
  1264.     /* read in parameters */
  1265.  
  1266.     a = vec (MAX_PARAMS);
  1267.     if ((par_name = (fixstr *) malloc ((MAX_PARAMS+1)*sizeof(fixstr))) == NULL)
  1268.     Eex ("Memory running low during fit")
  1269.     num_params = 0;
  1270.  
  1271.     if (!equals(c_token, "via"))
  1272.         int_error("Need via and either parameter list or file", c_token);
  1273.     if (isstring(++c_token)) {
  1274.         boolean fixed;
  1275.         double    tmp_par;
  1276.         char c, *s;
  1277.         char sstr[MAX_LINE_LEN];
  1278.         FILE *f;
  1279.  
  1280.         quote_str(sstr, c_token++, MAX_LINE_LEN);
  1281.  
  1282.         /* get parameters and values out of file and ignore fixed ones */
  1283.  
  1284.         fprintf (log_f, "take parameters from file: %s\n\n", sstr);
  1285.         if ( !(f = fopen (sstr, "r")) )
  1286.             Eex2 ("could not read parameter-file %s", sstr);
  1287.         while ( TRUE ) {
  1288.             if ( !fgets (s = sstr, sizeof(sstr), f) )        /* EOF found */
  1289.                 break;
  1290.             if ( (tmp = strstr (s, FIXED)) != NULL ) {          /* ignore fixed params */
  1291.                 *tmp = '\0';
  1292.                 fprintf (log_f, "FIXED:  %s\n", s);
  1293.                 fixed = TRUE;
  1294.             }
  1295.             else
  1296.                 fixed = FALSE;
  1297.             if ( (tmp = strchr (s, '#')) != NULL )
  1298.                 *tmp = '\0';
  1299.             if ( is_empty(s) )
  1300.                 continue;
  1301.             tmp = get_next_word (&s, &c);
  1302.             if ( !is_variable (tmp)  ||  strlen(tmp) > MAX_VARLEN ) {
  1303.                 fclose (f);
  1304.                 Eex ("syntax error in parameter file");
  1305.             }
  1306.             
  1307.             copy_max (par_name[num_params], tmp, sizeof(fixstr));
  1308.             /* next must be '=' */
  1309.             if ( c != '=' ) {
  1310.                 tmp = strchr (s, '=');
  1311.                 if ( tmp == NULL ) {
  1312.                     fclose (f);
  1313.                     Eex ("syntax error in parameter file");
  1314.                 }
  1315.                 s = tmp+1;
  1316.             }
  1317.             tmp = get_next_word (&s, &c);
  1318.             if ( sscanf (tmp, "%lf", &tmp_par) != 1 ) {
  1319.                 fclose (f);
  1320.                 Eex ("syntax error in parameter file");
  1321.             }
  1322.             
  1323.             /* make fixed params visible to GNUPLOT */
  1324.             if ( fixed ) {
  1325.                 struct value v;
  1326.                 Gcomplex (&v, tmp_par, 0.0);
  1327.                 setvar (par_name[num_params], v);   /* use parname as temp */
  1328.             }
  1329.             else
  1330.                 a[num_params++] = tmp_par;
  1331.             
  1332.             if ( (tmp = get_next_word (&s, &c)) != NULL ) {
  1333.                 fclose (f);
  1334.                 Eex ("syntax error in parameter file");
  1335.             }
  1336.         }
  1337.         fclose (f);
  1338.     } else { /* not a string after via */
  1339.         fprintf (log_f, "use actual parameter values\n\n");
  1340.         do {
  1341.             if (!isletter(c_token))
  1342.                 Eex ("no parameter specified");
  1343.             capture(par_name[num_params], c_token, c_token, sizeof(par_name[0]));
  1344.             /* create variable if it doesn't exist */
  1345.             a[num_params] = createdvar (par_name[num_params], INITIAL_VALUE);
  1346.             num_params++;
  1347.         }    while (equals(++c_token, ",") && ++c_token);
  1348.     }
  1349.     redim_vec (&a, num_params);
  1350.     par_name = (fixstr *) realloc (par_name, (num_params+1)*sizeof(fixstr));
  1351.  
  1352.     /* avoid parameters being equal to zero */
  1353.     for ( i=0 ; i<num_params ; i++ )
  1354.     if ( a[i] == 0 )
  1355.         a[i] = NEARLY_ZERO;
  1356.  
  1357.     regress (a);
  1358.  
  1359.     fclose (log_f);
  1360.     log_f = NULL;
  1361.     free (fit_x);
  1362.     free (fit_y);
  1363.     free (err_data);
  1364.     free (a);
  1365.     free (par_name);
  1366.     if ( func.at ) {
  1367.         free (func.at);                     /* release perm. action table */
  1368.         func.at = (struct at_type *) NULL;
  1369.     }
  1370.  
  1371. }
  1372.  
  1373. #if defined(ATARI) || defined(MTOS)
  1374. int kbhit()
  1375. {
  1376.     fd_set rfds;
  1377.     struct timeval timeout;
  1378.     
  1379.     timeout.tv_sec = 0;
  1380.     timeout.tv_usec = 0;
  1381.     rfds = (fd_set)(1L << fileno(stdin));
  1382.     return((select(0, &rfds, NULL, NULL, &timeout) > 0) ? 1 : 0);
  1383. }
  1384. #endif
  1385.